Plaquette Renormalization Scheme for Tensor Network States 
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We present a method for contracting a square-lattice tensor network in two dimensions, based on 
auxiliary tensors accomplishing successive truncations (renormalization) of 8-index tensors for 2x2 
plaquettes into 4-index tensors. The scheme is variational, and thus the tensors can be optimized 
by minimizing the energy. Test results for the quantum phase transition of the transverse-field Ising 
model confirm that even the smallest possible tensors (two values for each tensor index at each 
renormalization level) produce much better results than the simple product (mean-field) state. 
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Tensor network states (TNSs)"=»i» are emerging as a 
promising route toward unbiased modeling of challenging 
quantum many-body systems, such as frustrated spins. 
These correlated states are higher-dimensional general- 
izations of matrix product states (MPSs),— <2 which are 
implicitly produced in density matrix renormalization 
group (DMRG) calculations^ and are known to faith- 
fully represent ground states of one-dimensional (ID) 
hamiltonians with short-range interactions £ The matrix 
size to has to increase at most polynomially with the 
system size N, which underlies the success of the DMRG 
method in ID. For 2D and 3D systems, correlations are 
not reproduced properly by MPSs^ due to the inher- 
ently ID nature of the local quantum entanglement in 
these states (although improved scheme o 11 ' 12 can restore 
2D or 3D uniformity), and to then has to grow expo- 
nentially with N . In the TNSs, the matrices are re- 
placed by tensors of rank corresponding to the coordi- 
nation number of the lattice, e.g., on a 2D square lat- 
tice the tensors have four indices, in addition to their 
physical index (here spin z-component), as illustrated in 
Fig. Q] While it is believed, based on entanglement en- 
tropy considerations^ that TNSs can represent ground 
states of short-range 2D and 3D hamiltonians, a seri- 
ous problem in practice is that contracting the tensors 
is, in general, an NP hard problem. To overcome this 
challenge, approximate ways to compute the contraction 
have been proposed i 1 ' 13 ! 14 Another approach is to use 
tree-tensor networks^ or more sophisticated extensions 
of these,— which can be efficiently contracted. Promis- 
ing results based on TNSs have already been reported 
for several quantum spin models, but further reduction 
of the computational complexity, while maintaining the 
ability of the TNSs to properly account for entanglement, 
will still be necessary before the most challenging systems 
can be studied reliably. 

In this paper, we introduce a plaquette renormaliza- 
tion scheme for 2D TNSs inspired by recent work of Levin 
and NaveJ^ They suggested to replace the effective ten- 
sors for 2 x 2 plaquettes on the square lattice, which have 
m 8 elements, by some "renormalized" tensors with m^ ut 
elements, as illustrated in Fig. [DJa). This is exact, cor- 
responding just to a regrouping of indices, if m cut = to 2 , 



and may be a good approximation even if m cu t ~ to. For 
a square L x L lattice, the new tensors describing 2x2 
spins should be contracted on a new lattice of length L/2. 
If the original L is a power of 2, this decimation can be 
continued until there is a single tensor left, which is then 
contracted with itself (under periodic boundary condi- 
tions) to give the wave function. One can also carry out 
this kind of renormalization of the "double tensors" ob- 
tained when the physical spins are traced out first. The 
full contraction of the double-tensor network, illustrated 
in Figs.QJbjc), is the norm (\I r |\I r ) of the wave function. A 
matrix element (Xl'lyll 11 !/) of some operator involving one 
or several sites can be treated in a similar way. In prac- 
tice, in most TNS approaches, one would first construct 
the double tensors and then compute the full contrac- 
tion in some approximate way. A plaquette renormaliza- 
tion with double tensors is depicted in Fig. [2)Jb). This 
may be a good approximation for some judiciously cho- 
sen renormalized tensor, but how to find the best one 
is an open question. Gu, Levin, and Wen implemented 
a singular value decomposition (SVD) scheme^ for the 
double tensors, and a similar method was proposed by 
Jiang, Weng, and XiangiH In our scheme, the renormal- 
ization is instead accomplished with the aid of auxiliary 
3- index tensors S™ bc , which transform and truncate pairs 
of indices of the plaquette tensors in the wave function, 
as shown in Fig. (3J The sequence of plaquette renormal- 
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FIG. 1: (a) For each site s on an iV-site square lattice, there 
is a 4-index tensor T^ kl {a s ), where a s is the physical in- 
dex (here a spin a% = ±1) and i,j,k,l £ {l,...,m}. The 
wave function 9(<Ji, . . . , crjv) is the JV-tensor product (con- 
traction) defined as a summation over all shared indices on 
the links, (b) The norm |^) is obtained by contracting also 
the physical indices, which is normally done by first summing 
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FIG. 2: (a) Truncation (renormalization) of an 8-index pla- 
quette tensor into an effective 4-index tensor. The indices 
a, . . . ,h G {1, • • • , m} and i,j,k,l S {1, . . . , m cut }. The diag- 
onal lines indicate physical indices, which for an S — | spin 
system can take two values before the renormalization and 16 
values after (and is further multiplied by 16 after each succes- 
sive renormalization). In (b) the physical indices are traced 
out first, leading to double tensors. To stay with the same 
level of truncation as in (a), the indices after the renormal- 
ization should then take values i,j,k,l € {1, . . . , m 2 ut }. 

izations effectively corresponds to a new tensor network, 
which is illustrated in the case of the wave function on 
an 8 x 8 lattice in Fig. [4] 

Here we will apply this scheme to the transverse-field 
Ising model (which has become the bench-mark of choice 
for initial tests within TNS approaches); 

H = -J^M-hJ2<, CD 

(ij) « 

where erf and erf are standard Pauli matrices and (ij) 
denotes nearest-neighbor site pairs on a 2D square L x L 
lattice with periodic boundary conditions. The ground- 
state wave function of this system is translationally in- 
variant and positive-definite. We can then take the orig- 
inal tensors to be site-independent, i.e., there are just 
two tensors Tijki(<J s — ±1) and the tensors S n are the 
same for all plaquettes on a given renormalization level 
n = 1, log 2 (i) — 1. At the last level, four tensors 
remain (see Fig. [3]), which we contract directly. The 
problem is now to optimize the elements of the T and 
S tensors, to minimize the energy. Before proceeding to 
calculations, several comments are in order. 

It is clear from Fig. 0] that the plaquette renormaliza- 
tion (no matter how it is accomplished) breaks transla- 
tional symmetry, which, in the optimized state, should 
gradually be restored with increasing m. A way to re- 
store the symmetry for finite m cu t is to sum over all (here 
four) symmetrically non-equivalent ways of arranging the 
plaquettes on the lattice at each level. 12 In a similar way, 
one can also ensure that the wave function is symmet- 
ric under other lattice transformations (rotations and re- 
flections) for arbitrary T and S n (as an alternative to 
enforcing these symmetry in the individual tensors), and 
spin-inversion symmetry can be implemented in a similar 
way. These symmetrization procedures (which also en- 
able studies of states with different quantum-numbers of 



2 




FIG. 3: Renormalization of an 8-index plaquette tensor using 
auxiliary 3-index tensors. The solid circles denote either the 
original tensors T = T°, which depend on the physical spins 
(which are not shown here), or tensors T n arising after n 
renormalization steps have been carried out (and depend on 
the 2 2 ™ spins within the block). The squares denote 3-index 
tensors S n by which the external indices are decimated by 
contracting common T n_1 and S n indices. The remaining 
four free indices take values 1, . . . , m cut . 




FIG. 4: The effective reduced tensor network for an 8 x 8 
lattice. Here there are two sets of S tensors; S 1 and S 2 , 
denoted by the smaller and larger squares. The solid circles 
correspond to the original, spin dependent tensor T°. 

the symmetry operators) can be carried out if the spins 
are sampled using Monte Carlo simulations^ but can- 
not be easily used with the double-tensor network. Here, 
in this initial test of the scheme, we will trace out the 
spins exactly and work with the double tensors. When 
minimizing the energy, it is then important to calculate 
the full, translationally averaged energy, not just the site 
and bond energies on a single plaquette (as is done in the 
SDV scheme of Gu et alr^-), in order to maintain the vari- 
ational property of the scheme. No approximations are 
made when contracting the effective renormalized tensor 
network, also when using the double tensors. In contrast, 
the SVD applied to the double-tensor network^ intro- 
duces an approximation due to which the calculated en- 
ergy becomes non-variational. This can cause problems 
when optimizing the tensors. Note that in our double- 
tensor approach, there is a pair of equal S tensors for each 
plaquette edge (one from (\1/| and one from |<3>)) and one 
cannot combine these into arbitrary tensors, which would 
make the scheme non-variational. 

The renormalization of a plaquette according to Fig. [3] 
requires 12 internal index summations for each combi- 
nation of the four external indices of the renormalized 
tensor, which can be carried out with ~ m 8 operations 
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FIG. 5: Steps in the contraction of the renormalized plaquette 
tensor of Fig. [3] Partial contractions are constructed from left 
to right. The summations carried out are indicated by £ and 
the scaling of each step with m is shown beneath the diagrams. 
Each summation and free index contributes a factor m. 

(assuming m cut — m, which we will use here). One ex- 
ample of a sequence of contractions giving this scaling 
is shown in Fig. [5] When working with double tensors, 
the scaling becomes to 16 , since all external and internal 
indices should then take m 2 ut = to 2 values. The deriva- 
tives of the energy with respect to all the tensor elements, 
which we use to minimize the energy, can be evaluated 
with the same scaling in to, using a chain-rule procedure 
carried out along with the renormalization steps (as we 
will describe in detail elsewhere) 

It should be noted that the optimized T tensors do not 
necessarily constitute a good TNS when contracted with- 
out the S tensors, because they are optimized together. 
On the other hand, it should also be possible to construct 
special S tensors that effectively perform something very 
similar to a SVD (although globally optimized, not lo- 
cally as in and then the optimized T tensors 
by themselves should also form a good TNS when assem- 
bled into a standard 2D tensor network. Here we do not 
impose any such conditions on the S tensors. 

We now discuss calculations for the hamiltonian {!]). 
For L — > oo, the ground state of this system undergoes 
a quantum phase transition in the 3D Ising universality 
class, at a critical field h c /J w 3.04 determined using 
quantum Monte Carlo simulations^ Here we study the 
model using the smallest possible tensors and truncation, 
m = m cn t = 2, using the double-tensor approach for 
lattices of size L = 4,8, 16, and 32. For the T tensors, we 
enforce symmetry with respect to rotations of the indices, 
for a total of 6 free parameters each for T(l) and T(— 1). 
The S tensors are symmetric in the two indices used in 
the internal contractions of the plaquettes, and so there 
are 6 free parameters also for each S n . 

Starting with random tensor elements, we calculate 
the energy and its derivatives with respect to all the pa- 
rameters. We then update the parameters based on the 
derivatives, using a stochastic optimization scheme of the 
kind discussed in Ref.[l^. The plaquette renormalizations 
and the derivative calculations are highly parallelizable, 
which we take advantage of by using a massively parallel 



computer. In general the method performs very well, al- 
though occasionally the energy converges to a local mini- 
mum, especially close to the critical point. There can also 
be problems with the wave-function norm becoming very 
small, which can be remedied by periodically multiplying 
the tensors by suitable factors. 

The method can produce solutions breaking spin- 
inversion symmetry, in which case we do not obtain the 
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FIG. 6: (Color on-line) Spontaneous magnetization m z (up- 
per panel) and field-induced polarization m x (lower panel) 
versus the transverse field for different L x L lattices. 

true ground state for a given system size L, which can 
have no broken symmetries for finite L. This is to be 
expected, in analogy with how mean-field theory (corre- 
sponding to the to — 1 simple product state) produces 
symmetry-broken states for h < h c . Thus we can simply 
compute the magnetization m z = (<r|) (averaged over all 
sites s) and study its behavior for increasing L. Note that 
for fixed L, spin-inversion symmetry should gradually be 
restored with increasing m, as the optimized state should 
approach the true finite-!/ symmetric ground state. How- 
ever, for any fixed m, we expect the symmetry to be bro- 
ken when L — ► oo for h below some m-dependent h c . 
Here we only consider m = 2. Results for m z and the 
field-induced m x = (ctJ) are shown versus the field h/J 
in Fig. [5] The transition point between the magnetic and 
paramagnetic states moves toward higher fields with in- 
creasing L, converging to h c /J sa 3.33 for the largest 
sizes. This is closer to the unbiased quantum Monte 
Carlo result h c /J ~ 3.04 than the mean-field (to = 1) 
value h c / J = 4. There is some weak finite-size round- 
ing close to the transition for h > h c (L), which becomes 
less pronounced with increasing L. One the other hand, 
the magnetizaton curve for h < h c (L) becomes less sharp 
with increasing L. For L — 4 and 8 the transitions look 
almost first-order. However, upon closer examination, 
the curves can be described in a range of fields h < h c by 
a power-law, (h c — h)P, with the exponent (5 increasing 
with L. For L = 32 and 64, f3 = 0.40 describes the data 
well over a wide range of fields, but it is not clear whether 
this is the true asymptotic value. Including only points 
very close to h c , the mean-field exponent ft = 1/2 also 
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FIG. 7: (Color on-line) Spin-spin correlation versus sepa- 
ration r at fields h close to h c such that the long-distance 
correlations are approximately the same for system sizes 
L = 16, 32, and 64. 

works well. It is thus presently unclear whether our TNS 
can describe non-mean-field critical correlations for finite 
m = m cut = 2, although we expect that such behavior 
should certainly emerge with increasing m = m cut . 
Fig. [7] shows the spin-spin correlation function 



C(r 



{SfS*} averaged over all equidistant spins close 



to h c . This shows behavior clearly different from a sim- 
ple m = 1 product state, which only gives a constant 
C(r) — (Sf) 2 for r > 0. Characteristic boundary en- 
hancements of the correlations at r ~ L/2 are seen for 
the smaller sizes. 

In summary we have presented a scheme using auxil- 
iary tensors to renormalize plaquette tensors in a 2D ten- 
sor network. The approach can also be regarded as a dif- 
ferent tensor network, which can be contracted efficiently. 
Figs. [3] and 2] summarize the approach pictorially. Using 
the example of the transverse-field Ising model, we have 
shown that the scheme produces results far better than 
mean-field theory, even with the smallest possible non- 
trivial tensors and truncation (m = 2). The scaling to 
the thermodynamic limit is well-behaved. Based on these 



results, we expect a fast convergence to the exact ground 
state with increasing m. Increasing m to 3 is already 
quite challenging within the double-tensor approach that 
we have employed here, since the scaling is m 16 . How- 
ever, in variational Monte Carlo simulations (sampling 
the spins and optimizing the tensors based on stochas- 
tic estimates of the derivatives^) the scaling is m 8 , and 
it should then be easier to study larger m. Application 
to other quantum spin systems (and even fermions) is in 
principle straight-forward, although the convergence with 
m can of course be expected to be model dependent. 

Here we optimized the tensors variationally, which re- 
quires the energy averaged over all non-equivalent sites 
and bonds. The computational effort then scales with 
the system size as L 2 log(L). If we optimize only a local 
energy, which is not a variational estimate of the total en- 
ergy but can produce good results in the SVD approach 13 
(and may work well also in our scheme for some particu- 
lar classes of S tensors), the scaling is log(L), as in SVD 
based schemes. It may also be possible to use imaginary- 
time evolution (ground-state projection), as is often done 
in other TNS approaches iSS 

A technically appealing feature of our scheme is that 
the plaquette renormalization procedure can be imple- 
mented very efficiently on CPUs (graphics processing 
units, the use of which is emerging as an important 
trend in high-performance computingSi) . The speed-up 
relative to a standard CPU can be very significant for 
large tensors. For the spin-dependent tensors we have 
achieved better than a factor of 150 in efficiency boost 
when 7ti ~ 6. We plan to use this approach in future 
model studies. 
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